# Get Table 4
# =========================================================================================================
# remove everything
rm(list=ls())

# =================================================================================

# Load packages
library(foreign)

# Load the bootstrapped data
load("Boot_s_Dt.RData")
load("Jack_s_Dt.RData")

# ---------------------------------------------------------------------
# Get the bootstrapped CI.
B <- 400

# summary of the estimates
DLE <- summary(BpF[,1])
# summary of the Bootstrapped estimates
DLB <- apply(BpF, MARGIN = 2, summary)

# acceleration factor
DLE.jk <- apply(JKpF, MARGIN = 2, summary)
acc.top <- apply((apply(DLE.jk, MARGIN = 1, mean) - DLE.jk)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(DLE.jk, MARGIN = 1, mean) - DLE.jk)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# ---------------------------------------------------------------------
# 1st quartile: the 2nd position of DLE is the mean estimate
pos <- 2
dec <- 4
phio <- qnorm(sum(DLB[pos,] < DLE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(DLB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(DLB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q1 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Median: the 3rd position of DLE is the mean estimate
pos <- 3
dec <- 4
phio <- qnorm(sum(DLB[pos,] < DLE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(DLB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(DLB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.md <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# 3rd quartile: the 5th position of DLE is the mean estimate
pos <- 5
dec <- 4
phio <- qnorm(sum(DLB[pos,] < DLE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(DLB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(DLB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q3 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Percentage of positive DL obs.
dec <- 4
phio <- qnorm(apply((BpF[,1:400] < BpF[,1]), MARGIN = 1, sum)/B)

phi5 <- qnorm(0.05)

acc.top2 <- apply((apply(JKpF, MARGIN = 1, mean) - JKpF)^3, MARGIN = 1, sum)
acc.bot2 <- 6*(apply((apply(JKpF, MARGIN = 1, mean) - JKpF)^2, MARGIN = 1, sum))^(3/2)
acc.factor2 <- acc.top2/acc.bot2
c <- acc.factor2
a12 <- pnorm(phio + (phio+phi5)/(1-c*(phio+phi5)))
up <- c()
for(jj in 1:nrow(BpF)){
  temp <- as.numeric(quantile(BpF[jj,1:400], probs = a12[jj]))
  up <- rbind(up, temp)
  if(jj==round(jj,-3)) print(jj)
}
per <- format(round(sum(up>0)/nrow(BpF), dec), nsmall = dec)
print(per, quote = F)

# ---------------------------------------------------------------------
# Tabulate
print(paste(ci.q1, ci.md, ci.q3, sep = ";"),quote = F)

# ================================================================================================
# ================================================================================================

# summary of the estimates
SPE <- summary(Bpsw[,1])
# summary of the Bootstrapped estimates
SPB <- apply(Bpsw, MARGIN = 2, summary)

# acceleration factor
SPE.jk <- apply(JKpsw, MARGIN = 2, summary)
acc.top <- apply((apply(SPE.jk, MARGIN = 1, mean) - SPE.jk)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(SPE.jk, MARGIN = 1, mean) - SPE.jk)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# ---------------------------------------------------------------------
# 1st quartile: the 2nd position of SPE is the mean estimate
pos <- 2
dec <- 4
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q1 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Median: the 3rd position of SPE is the mean estimate
pos <- 3
dec <- 4
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.md <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# 3rd quartile: the 5th position of SPE is the mean estimate
pos <- 5
dec <- 4
phio <- qnorm(sum(SPB[pos,] < SPE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(SPB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(SPB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q3 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Percentage of positive SP obs.
dec <- 4
phio <- qnorm(apply((Bpsw[,1:400] < Bpsw[,1]), MARGIN = 1, sum)/B)

phi5 <- qnorm(0.05)

acc.top2 <- apply((apply(JKpsw, MARGIN = 1, mean) - JKpsw)^3, MARGIN = 1, sum)
acc.bot2 <- 6*(apply((apply(JKpsw, MARGIN = 1, mean) - JKpsw)^2, MARGIN = 1, sum))^(3/2)
acc.factor2 <- acc.top2/acc.bot2
c <- acc.factor2
a12 <- pnorm(phio + (phio+phi5)/(1-c*(phio+phi5)))
up <- c()
for(jj in 1:nrow(Bpsw)){
  temp <- as.numeric(quantile(Bpsw[jj,1:400], probs = a12[jj]))
  up <- rbind(up, temp)
  if(jj==round(jj,-3)) print(jj)
}
per <- format(round(sum(up>0)/nrow(Bpsw), dec), nsmall = dec)

# ---------------------------------------------------------------------
# Tabulate
print(paste(SPE[2],SPE[3],SPE[5],per, sep = ";"), quote = F)
print(paste(ci.q1, ci.md, ci.q3, sep = ";"),quote = F)
per

# ================================================================================================
# ================================================================================================

Bptil <- Bptil[complete.cases(Bptil),]; JKptil <- JKptil[complete.cases(JKptil),]
# summary of the estimates
TILE <- summary(Bptil[,1])
# summary of the Bootstrapped estimates
TILB <- apply(Bptil, MARGIN = 2, summary)

# acceleration factor
TILE.jk <- apply(JKptil, MARGIN = 2, summary)
acc.top <- apply((apply(TILE.jk, MARGIN = 1, mean) - TILE.jk)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(TILE.jk, MARGIN = 1, mean) - TILE.jk)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# ---------------------------------------------------------------------
# 1st quartile: the 2nd position of TILE is the mean estimate
pos <- 2
dec <- 4
phio <- qnorm(sum(TILB[pos,] < TILE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(TILB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(TILB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q1 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Median: the 3rd position of TILE is the mean estimate
pos <- 3
dec <- 4
phio <- qnorm(sum(TILB[pos,] < TILE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(TILB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(TILB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.md <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# 3rd quartile: the 5th position of TILE is the mean estimate
pos <- 5
dec <- 4
phio <- qnorm(sum(TILB[pos,] < TILE[pos])/B)
phia <- qnorm(0.025); phina <- qnorm(0.975)

# a_1 and a_2
c <- acc.factor[pos]
a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
up <- format(round(quantile(TILB[pos,], probs = c(a1)), dec), nsmall = dec)
down <- format(round(quantile(TILB[pos,], probs = c(a2)), dec), nsmall=dec)
ci.q3 <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")

# ---------------------------------------------------------------------
# Percentage of positive DL obs.
dec <- 4
phio <- qnorm(apply((Bptil[,1:400] < Bptil[,1]), MARGIN = 1, sum)/B)

phi5 <- qnorm(0.05)

acc.top2 <- apply((apply(JKptil, MARGIN = 1, mean) - JKptil)^3, MARGIN = 1, sum)
acc.bot2 <- 6*(apply((apply(JKptil, MARGIN = 1, mean) - JKptil)^2, MARGIN = 1, sum))^(3/2)
acc.factor2 <- acc.top2/acc.bot2
c <- acc.factor2
a12 <- pnorm(phio + (phio+phi5)/(1-c*(phio+phi5)))
up <- c()
for(jj in 1:nrow(Bptil)){
  temp <- as.numeric(quantile(Bptil[jj,1:400], probs = a12[jj]))
  up <- rbind(up, temp)
  if(jj==round(jj,-3)) print(jj)
}
per <- format(round(sum(up>0)/nrow(Bptil), dec), nsmall = dec)
print(per, quote = F)

# ---------------------------------------------------------------------
# Tabulate
print(paste(ci.q1, ci.md, ci.q3, sep = ";"),quote = F)




